%this script will find a non power positive Leslei matrix A
find=false;
while(~find)
    A=zeros(4);
    for n=1:4
        A(1,n)=fix(rand(1)*3);
        if(n<4)
            A(n+1,n)=roundn(rand(1)*0.9+0.1,-1);
        end;
    end;
    if(~IsPowerPositive(A))
        [V,D]=eig(A);
        if(imag(max(max(D)))~=0)
            find=true;
        end;
    end;
end;

%this script will find a power positive Leslei matrix B
find=false;
while(~find)
    B=zeros(4);
    for n=1:4
        B(1,n)=fix(rand(1)*3);
        if(n<4)
            B(n+1,n)=roundn(rand(1)*0.9+0.1,-1);
        end;
    end;
    if(IsPowerPositive(B))
        find=true;
    end;
end;

%display A and B
disp('This is the power positive matrix (B)');
disp(B);
disp('This is the non power positive matrix (A)')
disp(A);

%Calculate the minimum entry of A powers n=1:10
for n=1:10
    disp(strcat('The minimum entry of A^',int2str(n),' is:'));
    disp(min(min(A^n)));
end;

%Example of a non power positive growth
YearLimit=50;
Population=[[1;0;0;0],zeros(4,YearLimit)];
Popsum=zeros(1,YearLimit+1);
Popsum(1)=1;
for n=1:YearLimit
    Population(:,n+1)=A*Population(:,n);
    Popsum(n+1)=sum(Population(:,n+1));
end;

figure;
hold on;
plot(1:YearLimit+1,Population(1,:)./Popsum,'r');
plot(1:YearLimit+1,Population(2,:)./Popsum,'g');
plot(1:YearLimit+1,Population(3,:)./Popsum,'c');
plot(1:YearLimit+1,Population(4,:)./Popsum,'y');
xlabel('Time');
ylabel('Stage Percentage');
legend('Stage 1','Stage 2','Stage 3','Stage 4');